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ABSTRACT 

The uncertainty in the redshift distributions of galaxies has a significant potential impact on the 
cosmological parameter values inferred from multi-band imaging surveys. The accuracy of the pho- 
tometric redshifts measured in these surveys depends not only on the quality of the flux data, but 
also on a number of modeling assumptions that enter into both the training set and spectral energy 
distribution (SED) fitting methods of photometric redshift estimation. In this work we focus on the 
latter, considering two types of modeling uncertainties: uncertainties in the SED template set and 
uncertainties in the magnitude and type priors used in a Bayesian photometric redshift estimation 
method. We find that SED template selection effects dominate over magnitude prior errors. We 
introduce a method for parameterizing the resulting ignorance of the redshift distributions, and for 
propagating these uncertainties to uncertainties in cosmological parameters. 
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1. INTRODUCTION 

The surprisin g discovery of the accelerated expansion 
of the universe (jRiess et al.lll998HPerlmutter et al.lll999f ) 
has invigorated efforts to determine the history of the ex- 
pansion rate to high precision. By providing a greater 
understanding of the cosmic acceleration, further ob- 
servational input holds exciting prospects for improving 
our understanding of the fundamental laws of physics as 
well as the fate and possibly the origin of the universe 
(jAlbrecht et al.l 120061 ) . Many observational efforts are 
aimed at determining the distance-redshift relation via 
a variety of distance-dependent signals including the ap- 
parent magnitude of supernovae, the angular location of 
the acoustic feature in the galaxy number density corre- 
lation function, and the ampli tudes of tomographic cos - 
mic shear power sp ectra. See lBassett fc Hlozeld (2010 
Howell et al.l EMM : lleaTensl (|2009h : IZhan et al.l (pOtT 
Hutererl (|2010h for recent reviews. 

Determining the redshifts associated with these sig- 
nals can be challenging. In many cases, due to the large 
numbers and the faintness of the objects being used for 
distance determinations, we must extract redshift infor- 
mation not from spec troscopy, but from photometry i n 
multiple broad bands (Baum 1962; C onnolly et al.l l995). 
Redshifts determined from such data are called photo- 
metric redshifts, or simply photo-z's. Recent work has 
emphasized the stringent requirements on the quality of 
redshift information needed to avoid significantly degrad- 
ing or biasing inferences abo ut the dark energy posited to 
explain cosmic accelerat ion (Fernand ez-Soto et al.l 120021: 



Huterer fc Takadal l200l IMa et all 12006: 



Albrecht et al 



Bridle & King 



20061 : IZhanll20ul lAbdalla et al.l 12001 
20071 IMa fc Bernsteinll2008f ). 

Traditional tomographic analysis of photometric data 
for cosmic shear and galaxy clustering begins with a sort- 
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ing of galaxies into photo- z bins. This binning is a useful 
stage of data compression which, if done correctly, leads 
to very little loss of information. These bins can b e smal l 
in number and hence fairly coarse. As IMa et al.1 ff2006h 
have shown (in the context of a low-dimensional dark en- 
ergy model and a relatively shallow survey) cosmological 
parameter uncertainties do not decrease significantly in 
going from five to more photo- z bins. However, signifi- 
cant information loss can occur if the redshift probability 
distribution for a given galaxy, -P(z), is discarded imme- 
diately after binning. 

The analysis procedure we present here also begins 
with a sorting of galaxies. We refer to this sorting process 
as 'collecting' rather than 'photo- z binning' to emphasize 
two points: 1) the sorting criterion, (or possibly criteria) 
is not restricted to redshift estimates and 2) even if the 
sorting is based on a redshift estimate, the distribution of 
the collection of galaxies with redshift will extend beyond 
the photo- z redshift range used to define the collection. 

However the collection is defined, a weak lensing map 
and a galaxy count map can then be made for each col- 
lection. From these maps, all the auto and cross-power 
spectra can then be calculated (|Hulll999| ), as well as other 
summ ary statistics such as bispectra dTakada fc JaXnl 
|2004( ) and shear peak counts (jWang et al.ll2009f h 

The crucial redshift information that we need in or- 
der to calculate model predictions for these power spec- 
tra (and any other summary statistics) is the redshift 
distribution, dN/dz , for each co llection (|Ma et al.l 120061 : 
IHuterer et al.1 120061: IZh an 2006). In this paper, as dis- 
cussed in ICunha et alj (|2009f ) and IWittmanf (|2009D . we 
sum individual P(z)'s for the galaxies in each collection 
in order to infer the dN/dz for that collection. 

Our estimates of dN/dz are model- dependent. They 
depend on assumptions about the distribution of galax- 
ies with respect to spectral type, apparent magnitude in 
a reference band and redshift. We analytically demon- 
strate that, if the modeling assumptions are correct that 
go into our dN/dz estimates, then these estimates are 
unbiased (Appendix IA")) . 

Of course the modeling assumptions cannot be per- 
fectly correct. Thus a crucial element of our method- 
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ology is a treatment of the uncertainty in dN/dz that 
arises from uncertainty in the modeling. Our treatment 
allows for this uncertainty to be propagated through to 
the rest of the datwea reduction process. The method for 
determining the uncertainty in dN/dz in a way that can 
be fully propagated to scientific conclusions is the chief 
contribution of this paper. 

Past work on propagation of redshift errors through 
to the cosmological parameters has typically relied 
on highly idealized models for the distribution of the 
redshift errors. The redshift distributions considered 
are often far simpler than the distributions expected 
for real photometr ic surv eys (Ilbert, et al. 2009). For 
instance IMa et al.1 (J2006) ' assume the distribution of 
photo-z estimate for given true redshift, P (z p h\z trU e) , 
to be Gaussian, paramete rized by a mean and scatter 
as a function of redshift. iHuterer et al.l (|2006| ) assume 
that photo-z errors take the form of a bias that varies 
with z. In reality, degeneracies in color space for 
multiband photometric determination of rcdshifts can 
result in redshift error distributions that are asym- 
metrical and multi-peaked, with peaks separated by 
Az > 1. As demonstrated by a growing body of 
research, it is important to take into account these 
catastrophic photo- z errors and to add complexity 
to th e photo- z uncertaint y mode ls dSchneider et al. 
2006t lAmara fc Refregierl 120081: IMa fc Bernstein 



2008t ISunetaLl 120091 iBernstein fc Hutered 120101: 
Nishizawa et al.ll2010HHearin et alJl2010f ). 

Although sufficiently robust spectroscopic sampling 
may reduce the need for modeling assumptions and the 
resulting uncertainties, for proposed large-scale future 
weak lensing surveys such as LSST, PanSTARRS and 
DES, there is much to be gained from using very large 
numbers of faint and distant galaxies for the distance de- 
terminations. In the case of LSST, several billion galaxies 
with i band magnitudes as faint as 25 will be used. It is 
very challenging to spectroscopically determine redshifts 
for a fair sample of such faint objects. Spectroscopic 
samples preferentially contain galaxies with clear spec- 
tral features and are therefore biased towards specific 
galaxy types. Moreover they are significantly volume 
incomplete several magnitudes brighter than the faint 
limits of the photometry In other words, training sets 
used to develop photometric model inputs may be un- 
representative of the larger galaxy population, resulting 
in uncontrolled errors in the photo-z estimation. Thus, a 
method for quantifying such modeling uncertainties is an 
absolute necessity for understanding the resulting uncer- 
tainties in the redshift distributions and the subsequent 
uncertainties in the cosmological parameters. 

We present here a method for parameterizing the red- 
shift error distributions that is completely general and 
thus accommodates the existence of catastrophic errors. 
The parameterization arises naturally out of the consid- 
eration of specific modeling assumptions and data. For 
specificity we use the LSST data model (described in §[3]) 
when giving examples of sensitivity to model assump- 
tions. 

We show that changes in modeling assumptions result 
in variations in the redshift distribution, dN/dz, which 
can be captured by principal component analysis (PC A), 
where just a few modes are sufficient to represent these 
effects. The dN/dz can be reproduced from linear com- 



binations of these modes with the weights on each mode 
becoming the parameters of the new representation. The 
uncertainty on each weight is extracted from the varia- 
tions resulting from uncertainties in modeling assump- 
tions, and this information is then usable for any like- 
lihood analysis into which photo-z uncertainties enter. 
These new variables are continuous and lend themselves 
to a Fisher matrix or MCMC analysis of the parameter 
space. 

In §[5] we review Bayesian (model-based) methods 
for inferring redshifts from multi-band photometry, and 
we introduce principal component analysis (PCA) as a 
method for naturally capturing the modeling uncertain- 
ties in photometric redshift estimation. In §[3] we present 
a model demonstrating how uncertainty in the selection 
of a set of spectral energy distribution (SED) templates 
can be captured by PCA. In §3] we explore uncertain- 
ties in the type and magnitude distributions tha t can be 
used as priors on photometric redshift estimates ()Benitezl 
2000) and demonstrate again the utility of PCA in pa- 
rameterization of the modeling uncertainty. In §[5] we 
relate our results to work that has propagated photo- 
z uncertainties through to constraints on dark energy. 
In §[5] we summarize our findings. 

2. INFERRING REDSHIFTS FROM MULTI-BAND 
PHOTOMETRY 

Our method takes as a starting point a calculation of 
the probability distribution of the redshift of an object 
given its measured colors and certain modeling assump- 
tions. We then present a simple estimator for the redshift 
distributions, dN/dz, of collections of such galaxies. For 
us, the "collections" will consist of galaxies all in the 
same photometric redshift bin. Of course, our model- 
ing assumptions may be incorrect. Thus we present a 
method for parameterizing the resulting uncertainty in 
the dN/dz for each photometric redshift bin. With the 
uncertainty in the dN/dz parameterized, and the uncer- 
tainty in the parameters specified, the problem is now 
well-posed for further propagation of the uncertainty in 
the photo-z 's to uncertainty in the cosmological param- 
eters. The following subsections lay out each of these 
steps. 

2.1. Bayesian Calculation of Redshift Probability 
Distributions 

Let us consider t to be a variable (or set of variables) 
that uniquely specifies the SED of a galaxy. Then spec- 
ifying for a galaxy its t, redshift z, and apparent mag- 
nitude in a single reference wavelength band, mi, is suf- 
ficient for determining the fluxes, fx, in a set of well- 
defined bandpasses where A = 1 we will take to be the ref- 
erence band. If we model the measured fluxes (the data), 
as due to the true fluxes plus noise terms f\=f\ + n\ 
and we assume that the n\ are Gaussian random vari- 
ables with mean zero and (n\n\> 
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We can turn this around, using Bayes' theorem, to 
get an expression for the probability distribution of the 
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unknowns t, z and mi in terms of the knowns, the ft, and 
whatever assumptions we care to make, which we shall 
call I (for Information). According to Bayes' theorem, 
which follows from axiomatic properties of probability 
distributions, 



P{t,z,m x \ftl) = P{fi\t,z, mi ,al) 



2 P(t,2,mi|J) 



P{ft\I) 



(2) 



The probability distribution on the left-hand side is 
called the "posterior" distribution (because it is the prob- 
ability distribution after we have the data), and the 
P{t, z,mi\I) is called the "prior" distribution (because 
it is the probability distribution for the parameters prior 
to our collection of the data). The prior is where all 
our modeling assumptions are encapsulated. The term 
in the denominator we can ignore since it has no depen- 
dence on the model parameters. While mi is formally a 
measured quantity derived from one of the band fluxes, 
it is treated as an unknown here because we want to map 
variations in P{z) by perturbing it within its Gaussian 
noise distribution. 

If we are interested in only the redshift, we can obtain 
its (marginal) distribution by integrating over all possible 
values of t and mi so that 



p{z) = p{z\f d x j) 



dt J dmxP(t,z,mi\ff,I) 



(3) 



where we have written the distribution simply as P(z), 
suppressing the dependence on the data and / (model- 
ing assumptions) for notational convenience^ Although 
we have ignored the P(f^\I) factor, we can now make 
urj^for that by choosing a normalization factor such that 
L P(z)dz = 1. A number of authors have used Bayes' 
theorem in this way to obtain the marginalized p robabil- 
ity d i stribution of the redshif t (see for instance iBenited 
[2000t lEdmondson et al.ll200& iChapin et all 12001 7" 

2.2. Estimating dN/dz 

Often the redshift distribution is reduced to a single 
value of redshift, z pi perhaps the highest peak in P{z), 
and then this value of z is referred to as a Bayesian photo- 
metric redshift estimate. Here we use z p to sort galaxies 
into photometric redshift bins (collections) but use all 
the information in the P(z) for each galaxy to estimate 
dN/dz for the collection. 

As various authors have noted (iBordoloi et all 120101: 
Cunha et all [2009; Mandelbaum et alJ l2008t iWittmanl 



20091 ). the use of the full P{z) can improve on the sys- 



tematic errors in photometric redshift estimation. Errors 
in photometric redshifts are often due to degeneracies in 
color space. (For instance, when an object's spectrum is 
represented by the fluxes in low resolution optical bands, 
the Lyman break at z ~ 3 may be indistinguishable from 
the 4000A break at low redshift.) Although this might 
result in a drastically miscalculated point estimate red- 
shift (i.e. catastrophic error), the P(z) for most galaxy 
types near this redshift is multimodal. Because inference 
of cosmological parameters (for instance through weak 
lensing observations) relies less on the redshift of specific 
galaxies and more on the distribution of galaxies in a 

4 Recall that t may be representing an array of variables and so 
the integration over t may be a multidimensional one. 
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Fig. 1. — The blue line is an example of a photometrically calcu- 
lated probability distribution for a single galaxy. The true redshift 
(the red 'x') is at z = 2.6, whereas the P(z) peaks at .2. By assign- 
ing non- negligible probability at the true redshift, it can be seen 
how the P(z) can be more accurate than a single-point estimate. 

collection, the use of the full P(z) is a vast improvement 
over point estimates. Consider Fig. [TJ the maximum 
of the P(z) is at low redshift, whereas the true redshift 
is near z ~ 3. In a point estimate approach this would 
be a catastrophic error, but the full P(z) assigns almost 
half of the probability to the galaxy being near its true 
redshift. 

Enumerating the galaxies in a given collection (redshift 
bin) with index g, we define our estimate for dN/dz in 
that collection as 



dN_ 

dz 



(z)="£P g (z). 



(4) 



The estimate is motivated by the fact that the probability 
distribution of the redshift of a galaxy drawn at random 
from the collection is given by jjTlgPgiz) where N is 
the number of galaxies in the collection. The AppendixlAl 
proves that this intuitive estimator is not biased as long 
as the prior and the likelihood in Equation (J - } arc known 
perfectly. 

The advantage to calculating the dN/dz in this way 
can be seen in Fig. [21 which shows the dN/ dz calculated 
for a mock galaxy catalog of 100, 000 galaxies consisting 
of the 6 CWWSB templates which follow the t ype dis- 
tribut ions as a function of magnitude given by iBenitezi 
(2000) (but with a fixed dN/dz peaking at z=2.0 for sim- 
plicity). The full P{z) is calculated for each galaxy, and 
the galaxies are sorted into redshift bins, or "collections" , 
based on the peak value of the P(z) (in other words, the 
photo-z point estimate). We compare the actual distri- 
bution, dN/dztme, to the dN/dz calculated from sum- 
ming and histogramming the single point photo-z esti- 
mates, dN/dzpczk- Because the calculated distribution, 
dN/dz pca k, is constructed from the same photo-z point 
estimates that were used to assign each object to a col- 
lection (bin) , the estimated redshift value necessarily lies 
within the redshift range for which the collection (bin) 
is defined; thus the calculated dN/dz pea k will always be 
zero outside that redshift range (or photometric "red- 
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Fig. 2. — dN/dz for two redshift "collections". Galaxies are 
sorted into collections on the basis of peak redshift value (or point 
estimate for the photo-z). The top plot shows collection 1, bottom 
plot shows collection 2. The blue line shows the dN/dz estimated 
by sorting galaxies into collections based on their most probable 
redshift and summing the single point photo-z estimates (the peak 
values of the P(z)). The circles are the value of the correspond- 
ing dN/dz for the true redshifts. The dashed black line shows the 
dN/ dz estimated by collecting galaxies using the peak of the prob- 
ability distribution and summing the entire P(z). It can be seen 
that dN/dz pro i, traces much more closely the true distribution than 
the distribution derived from the point estimate. 

shift bin"). This is a major failing of the point estimator 
method since the true redshift distribution of galaxies, 
dN/dztrue, may have tails that extend to redshifts above 
and below the redshift values for which the collection is 
defined. By using the full P(z), however, a galaxy placed 
in a collection can have a non-zero probability extending 
to values outside of redshift range of the collection. Thus 
the distribution calculated using the P(z) of each galaxy, 
dN/dZp TO b, can more accurately reflect the true redshift 
distribution. 

Despite the use of the full P(z) there are many possi- 
ble sources of error in the photo-z estimation that could 
still cause us to miscalculate the distribution of galaxies, 
dN/dz. These include use of a small number of galax- 
ies, mischaracterization of the statistical properties of 
the flux noise, and errors in the modeling assumptions 
summarized by the prior, P(t, z,mi\I). The latter in- 
cludes prior assumptions about galaxy distributions as 
a function of type, redshift, magnitude, surface bright- 
ness, etc., and also includes the choice of SED templates. 
Errors from these modeling assumptions are the central 
concern of this paper. 



2.3. Parameterizing Modeling Uncertainty 

We capture the effect of modeling uncertainty by ex- 
plicitly varying the modeling assumptions and seeing how 
the dN/dz's vary in response. We use PCA to define 
functions of z (modes) that describe the departures from 
the mean dN/dz. These modes, together with their am- 
plitudes, define our parameterization of dN/dz. 

We consider two different kinds of variation. In one 
procedure we explore effects due to discrete variations 
in the SED templates used to generate the photo-z es- 
timations. The other procedure is for the case where 
parameters governing the prior are not perfectly known, 
but are constrained by auxiliary observational data. In 
this latter case one can sample from the posterior dis- 
tribution of these "prior" parameters, and generate the 
dN/dz from this sample. 

No matter how one generates the samples of dN/dz, 
the PCA is performed as follows. Given a data matrix 
comprised of multiple entries for a set of observables, 
PCA analysis finds the eigenvalues and eigenvectors of 
the covariance matrix for these data. In this framework 
we treat the values of dN/dz at discrete points in red- 
shift space as our set of observables, and each modeling 
variation results in small changes in this set. We start 
by forming the data matrix 77 for each collection 



dN 

dz 



(5) 



Here i indexes the points in z space for which dN/dz is 
defined and a runs over the variations. Subtracting the 
mean of each row from each element in the row of 77, we 
have 



D 



where 



ia — r\icL r]{ 



Vi = —Y] Via 



(6) 



(7) 



and m is the number of variations. The covariance ma- 
trix is 

C = — - — DD T . (8) 

m — 1 

Computing the eigenvalues ({A} and eigenvectors {U}, 
the first principal component U\ corresponds to the 
largest eigenvalue, the second principal component U2 
corresponds to the next largest eigenvalue, and so on. 
Each data vector (indexed by a) can be written as a 
linear combination of the principal components, i.e. 



Vi 



N 

E 

13=1 



B a pUi/3 



(9) 



where the weights, B a p, are given by the linear transfor- 
mation 



B = U T r? 



(10) 



There are N principal components of length N, where 
again, N is the number of discretized values of z. If the 
eigenvalues decrease sufficiently quickly from biggest to 
smallest, then most of the variation in rj ia is described by 
the first few modes. By keeping only the first k modes, 
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we are able to reduce the dimensionality of the parameter 
space. We can reproduce the data via 

k 

vL = % + J2 B ^ u ^ ( n ) 

13=1 

where rj' ia ~ r\i a for a small value of k if the eigenvalue 
spectrum is sufficiently steep. 

The above description was for a single collection 
("photo-z bin") defined for a given redshift range. It 
can be repeated for each of the 7=1, ...n z photo-z col- 
lections; we will distinguish quantities for collections de- 
fined for different redshift ranges with a superscript 7. 

We can now finally define our parameterization of the 
redshift distribution as the amplitudes of the modes TJ7„ 
defined by the above process. We will call these am- 
plitudes b~p and they influence the redshift distributions 

rjj = dN/dz? via 

k 
/3=1 

Each modeling variation, a, leads to a particular value 
of bj — B^g and we can use these 'realizations' to 

model the prior distribution of the bZ, in this lower, k- 
dimensional space. The bj are not fixed at the discrete 
values, but are allowed to vary continuously to describe 
the full space of model uncertainty. This distribution 
will have zero mean due to the mean removal that occurs 
prior to the principal component decomposition. We can 
estimate a covariance matrix via 

1 m 

C^p, EE (b}b%) - I E V^'- (13) 

a=l 

2.4. Relation to the Estimation of Cosmological 
Parameters 

For the analysis of cosmic shear data, we assume it 
is possible to take a set of model shear power spectra, 
C^ v and calculate its likelihood, given the shear data. 
We further assume that one can take a specification of 
cosmological parameters, and a redshift distribution for 
each photo-z collection, and use these to calculate the 
model Cj 7 ". From the power spectra likelihood, any prior 
constraints on the cosmological parameters (for example 
from cosmic microwave background observations), and a 
prior constraint on the redshift distribution parameters, 
one can then form a joint posterior distribution for the 
cosmological plus redshift-distribution parameters: 

lnPpostcrior(^cos) = hi £ (shear datalC^^cos, bj)j 

+ lnP p rior (^cos )+lnP pri or(^) (14) 

where 

mP P ri OT (^) = Mw/J, W2 (15) 

With the ability to calculate this joint posterior, one 
would then be able to explore the constraints on this 
parameter space given by the data via Markov Chain 
Monte Carlo, for example. 



3. EXAMPLE 1: TEMPLATE UNCERTAINTY 

One of the potentially largest sources of modeling er- 
rors in photometric redshift estimation may come from 
the selection of the template set in the SED fitting 
method, the technique currently used for deep and faint 
redshift surveys. At present, many photo-z codes use 
a small number of template spectra fro m nearby galax- 
ies e . g. the popular C WW+SB set (jColeman et al.1 
1980; Kin ney et al.l 119961 ) used in part of our analysis. 
Though the SED s are sometimes 'twea ked' to better fit 
the data, (e.g. iBudavari et al.l [20001 ). expanding be- 
yond such simple template sets almost always involves 
the use of population s ynthesis models (|Brammer et al.l 
2008; Ilbe rt et al.ll2009l ). due to the lack of sizeable num- 
bers of high signal-to-noise spectra at high redshift. The 
difficulty of spectroscopic followup of faint, high redshift 
galaxies means these models may be incomplete and un- 
representative of true galaxy evolution. Such biases can 
have a signi ficant effect on the photo- z es timations (see 
for instance iMacDonald fc Bernsteinll2010h . Because the 
SED shape for some population of galaxies may evolve, 
or a population not represented at low redshifts/high lu- 
minosities is present in the deep photometric data, the 
choice of template set is an important source of modeling 
uncertainty. 

There are also convergence issues. Since the photo-z is 
a nonlinear inverse problem, a small amount of photom- 
etry noise will drive a noisy assignment of templates for 
each galaxy. The amount of this induced template noise 
is a nonlinear function of the photometry noise. This 
is a general feature of threshold systems: the signal-to- 
noise ratio peaks at a sma ll non-zero value of input noise 
(jWiesenfeld fc M oss 19951 ). At zero noise the photo-z er- 
ror is large and due entirely to errors in the modeling 
assumptions, whereas at high noise levels the photo-z er- 
ror is large because no template assignment converges. 
In this paper we focus on photo-z errors due to errors in 
models of SED templates and magnitude and type pri- 
ors at constant photometry noise given by the LSST data 
model. 

To begin to explore the sorts of uncertainties present 
in an actual survey with noise and where galaxy types 
may be left out or incorrectly modeled, we start by gen- 
erating a mock catalog based on 20 SED templates. We 
then calculate the photometric redshifts while system- 
atically leaving out one each of the SED's at a time. 
When one of the templates is removed, the galaxies in 
the catalog matching this SED will be fit to another, 
incorrect, SED, creating errors in the photometric red- 
shifts and variations in the resulting dN/dz. The result 
is a set of 21 distinct photometric redshift estimates for 
the galaxies in the catalog, which we can use to develop 
a parameterization and model of uncertainty for the red- 
shift distribution. (We note that this is just one of the 
many ways in which template noise can be examined.) 

3.1. Mock Data 

To generate our set of template s, we begin with acat- 
alog of galaxies from GOODSN (|Giavalisco et all 120041) 
with measured spectroscopic redshifts and seventeen 
band photometry, covering U-band to IRAC mid-infrared 
wavelengths, available as part of the Photo-z Accu- 



6 



racy and Testing (PHAT) program (jHildebrandt et al.l 
2010). Although there are hundreds of galaxies in the 
sample, for tractability we develop a more manageable 
set of 20 SED templates that are representative of the 
larger dataset. 

We derive o u r tem p late s et using the methodology 
of lAssef et al.l pOlOl . [20081) in which the 17 band 
data for each galaxy is fit with a non-negative linear 
combination of four basis templates (Elliptical, Spi- 
ral, Irregular, and active galactic nuclei (AGN)) shifted 
to the measured redshift. Modifying the code available at 
|http://www.astronomy.ohio-state.edu/~rjassef/lrt/index. 
to allow for the extraction of the coefficient associated 
with each component of each component and eliminating 
galaxies with a strong AGN component, we arrive at 
a sample of 375 galaxies, each described by a linear 
combination of E, Sbc, and Im templates. Examination 
of this reduced three dimensional space reveals that 
the galaxies' coefficients are not significantly clustered. 
In order to choose a more manageable, but still repre- 
sentative, set of templates that s pans the dataset, we 
run a simple k-means algorithm (Mac Queenl 119671 ) to 
choose 20 sets of coefficients that will comprise our mock 
galaxy template set. The resulting templates are shown 
in Fig. |31 To facilitate the following step, we categorize 
the templates in terms of their dominant morphological 
type, E, or Ti, Sbc or T 2 and Im, T3, with the reddest 
templates 1 — 6 corresponding to Ti, templates 7—13 
to T2 and the bluest templates 14 — 20 corresponding to 
T 3 . 

We generate magnitudes for the LSST ugrizy filter set, 
with depth and error properties expected for the "gold" 
sample described in section 3 of the LSST Science Book 
(LSST Sci ence Collaborations et all I2009D . which con- 
sists of expected full ten year depth and cut to i < 25.0 
magnitude. The 5a limiting magnitudes in the ugrizy 
bands are: 26.1, 27.4, 27.5, 26.8, 26.0, 24.8. The "gold" 
sample galaxies will have S/N > 20 in multiple filters, 
which is highly recommended for meaningful photo-z 
measurements. We begin by generating an r-band ap- 
parent magnitude that follows p(m) cx 10°- 37m to ap- 
proximate a typical number count distribution. We then 
choose one of the 20 template types for each galaxy by 
Monte Carlo, such that the morphological types, Ti, T 2 
and T3 are distributed according to 



TABLE 1 

Fiducial parameter values for magnitude priors 



Spectral type 


ft 


at 


zot 




kt 


E/SO 


0.35 


2.46 


0.431 


0.0913 


0.450 


Sbc,Scd 


0.50 


1.81 


0.390 


0.0636 


0.147 


Irr 


0.15 


0.91 


0.0626 


0.123 





html 



p(T\m ) = f t exp{-k t (m - 20)}. 



(16) 



We then assign via Monte Carlo a redshift to each galaxy, 
but such that 



is 1000 
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Fig. 3.— The 20 LRT SED templates derived from GOODSN 
and linear combinations of the Assef basis templates. The galaxies 
are further broken into three type classifications which follow the 
redshift distributions as parameterized in Table [T] 



p(z\T, mo) oc z a * exp{— [ 



(17) 




is satisfied, where z mt (mo) = zo t + fc mt ( m Q — 20), and 
T refers to E, Sbc, or Im. The free parameters, ft, kt, 
zg t and at, arc fit to a training set of spectroscopically 
m easured r edshif ts. We take these to be the values given 
in lBenitezI (|2000l ). (listed in Tabled]). Given the redshift, 
type, and r-band magnitude, the expected magnitudes in 
the remaining ugizy bands are assigned appropriately. 

5 Further information on PHAT, including the dataset used here 
is available at: http://www.astro.caltech.edu/twiki_phat/bin/view 
/Main/ WebHome 



Fig. 4. — LRT template sets: top panels show the dN/dz's from 
summing the P(z) of galaxies in each collection defined for a given 
redshift range. The bottom panels show the residuals for each 
realization - that is the difference between the fiducial (all 20) 
template set and the reduced (19 only of the 20) template sets, 
AdN/dzi = dN/dzfid — dN/dzi, for each template set realization. 

3.2. Results 

With the fiducial catalog generated, we can begin to 
explore how template noise affects the photo-z distribu- 
tion. We run the photo-z estimation code BPZ ([Benitezl 
2000) to calculate the P(z) for each galaxy in the cat- 
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alog. We perform 21 iterations - once with the fiducial 
model, e.g. all 20 SED's included for the estimation, and 
then 20 variations, each realized by leaving out one of 
the 20 SED templates. 

Each incomplete set of SED's produces slightly differ- 
ent errors in the photo-z estimations, which result in 
variations in the dN/dz for each collection. The dN/dz 
is stored as a vector of 134 "microbins" with a width of 
.03 from z = to z = 4. Hence the sum, 



134 



dN/d; 



'311 



(18) 



returns the total number of galaxies per collection 7, 
and summing over the galaxies in each collection re- 
turns the total number of galaxies in the catalog, e.g. 
100,000 = £ 7 iV 7 . In Fi g- H the dN/dz are shown for 
each realization of the templates. The bottom panels 
show the residuals, dN/dzi — dN/dz a u, where i refers to 
a realization where a single template has been left out for 
the photo-z estimation. It can be seen that the residuals 
are on the order of a few percent when compared to the 
dN/dz. 

The distribution of galaxies can be described by the 
mean and variance of the dN/dzi per collection 7, and 
are different for each template set realization. (Note that 
the collection index 7 has been suppressed for clarity.) 



(z)i 



and 



A,- 



dz z dN/dz (19) 



dz(z- (z)f 'dN/dz (20) 



where i refers to the specific template choice. We calcu- 
late the mean and the standard deviation of the previous 
quantities, i.e. 

1 ^ 

(21) 



» x = 21 2J Xi ) 



ox 



i 



1 21 



Xf 



(22) 



where X refers to either the mean or the variance of each 
variation. These values are shown in Fig. [5] (To keep 
the axes on the same scale, we actually plot the average 
offset for the bin, [i <z> - z cente r-) 

The spread in the means is typically on the order of 
10~ 3 , although for collection 4 it is larger. This reflects 
the difficulties in calculating photometric redshifts for the 
redshift range for collection 4 is defined, 2.4 < z < 3.2 
Specifically, the Lyman break at z ~ 3 is easily confused 
with the the 4000A at low redshift. (This is the cause 
of the low redshift "bump" in the dN/dz in collection 
4 seen in Fig. [4]) This effect is even more apparent in 
the standard deviation of the variance for collection 4. 
The variance is sensitive to probability moving into and 
out of outlying redshift "islands". To a lesser extent, 
this effect is also visible in the spread of the variance 
for collection 1. The effect is smaller, however, due to 
the large numbers of bright galaxies that also populate 
collection 1 and are never confused with distant galaxies. 
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Fig. 5. — LRT set: Top panel shows the average of the bias, 
/i< z > — z cen ter, where z C enter is the center of the redshift range 
for the collection. (This is essentially a plot of the average of 
the mean redshift for each variation, but the subtraction of the 
collection center allows the scale of top and bottom plots to be 
the same.) The standard deviation (spread) in the bias is plotted 
as error bars, where for visibility we multiply by a factor of 100: 
a X 100. The bottom panel shows the average variance, fj.^, and 
the standard deviation in the variance is plotted as error bars. For 
visibility we multiply by a factor of 10: a X 10. The large error bars 
in collections 4 and 5 show the effects of the degeneracies present 
in faint galaxies at low redshift. 

(Note the asymmetry of this degeneracy, which makes 
this effect more pronounced in the higher redshift bins: 
bright objects are unmistakably at low redshift, whereas 
faint objects may either be low redshift intrinsically faint 
galaxies or high redshift intrinsically bright galaxies.) 

We calculate the principal components of the dN/dz 
due to the variations of the template sets. The fraction 
of the variance captured by a given mode i is A Q J2i > 
where Xi is an eigenvalue of the covariance matrix corre- 
sponding to the ith mode. It is difficult to assign a phys- 
ical interpretation to the modes themselves, which show 
a complex and data-dependent structure. The modes 
will be a function of the methodology used to vary in- 
put model, e.g. we could have chosen to leave out 2 
templates at a time and the modes returned would have 
shown a different structure; however, they will always 
represent the "directions" in the data, rank-ordered by 
the amount of variance captured. As an example, the 
first three modes of collection 4 are plotted in Fig. [51 
Some features of the underlying data can be seen in the 
modes. The high amplitude of the modes at low redshift, 
for instance, reflects the fact that the template set vari- 
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Fig. 6. — LRT set: Top panel shows the principal components 
for collection 4 (illustrative of principal components for other col- 
lections, here simply labeled "bins"). The bottom panels show the 
eigenvalue spectrum (i.e. amount of variance captured by each 
principal component) (x-axis is logarithmic). In all cases, 80% or 
more of the variance is capture by the first three modes. 
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Fig. 7. — LRT set: Depicts d a u: the mean-subtracted dN/dz for 
the case of all 20 SED templates. The blue (solid) line is the the ac- 
tual residual. The cyan (dashed) line is the residual reconstructed 
from the first three principal components. Although some of the 
details (wiggles) are missed, the general shape is reproduced. 

ations move probability of certain galaxies into and out 
of this low redshift outlying region. 
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Fig. 8. — LRT set: The top panel shows the standard deviation 
of the variations d a (the mean-subracted dN/dz a ) at each point in 
redshift space. The bottom panel is the standard deviations of the 
d a reproduced from linear combinations of modes 4 and higher. It 
can be seen that most of the variance can be accounted for with 
modes 1 — 3 — little is described by the higher number modes. 

As noted above, the principal components are useful 
for dimensionality reduction. For all collections, the first 
three principal components cumulatively account for 80% 
or more of the variance. The eigenvalue spectrum is plot- 
ted in Fig. [6j where it can be seen that higher number 
modes quickly approach zero. 

Recall from Subsec. 12. 31 that the first step of PCA is to 
subtract the mean from each dN/dz a , where a refers to a 
specific realization of template set. We define {d a } to be 
the set of mean-subtracted data, i.e. dN / dz a — (dN / dz) . 
Each d a can be reconstructed from a specific linear com- 
bination of principal components. Performing the recon- 
struction of the dN/dz's of the LRT template variations 
with the first three modes, it is possible to recover 80% 
or more of the variance per collection. In general, the 
precise amount of variance recovered with the first three 
modes will depend on the specific methodology used to 
vary the photo- z modeling assumptions and also on the 
data set. In general, the variance and the fidelity of the 
reconstruction can be "tuned" by retaining a larger or 
smaller number of the modes, e.g. if 1% accuracy is re- 
quired, we can retain seven or eight modes that capture 
99% of the variance. In the case of three modes retained 
and five tomographic bins, the photo-z uncertainty is 
parameterized by just 15 free parameters. This is a rel- 
atively tractable number considering analyses propagat- 
ing redshift uncertainties through to the cosmol ogical pa- 
rameters have used 60 or more free param eters (|Ma et al.1 
l2?MIZr^l2?iMlMa fc Bernsteinll2008l) . 

A potential concern is that there is a mode (or modes) 
of dN/dz, discarded by the truncation, which, though not 
contributing much to variance of dN/dz, causes a change 
to the shear power spectrum (or other observable) with 
the same shape as caused by a cosmological parameter 
variation. Thus, although the amplitude of the variation 
is small, it could lead to bias in parameter estimation. 
We think such a scenario is unlikely, but one could test for 
it when this work is extended all the way to cosmological 
parameter estimation. 
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Additionally we note that knowledge of the specific 
redshift range of dN/dz to which a cosmological param- 
eter is most sensitive could improve the mode selection. 
Our PCA modes are defined and ordered based on the 
amount of variance they account for in dN/ dz. Of course 
the variation in the dN/dz is a function of redshift. For 
our analysis, we weight variations happening at differ- 
ent redshifts equally; however, one might improve the 
mode decomposition and their ordering by using some 
redshift-dependent weighting of the variation. We ex- 
pect the optimal weighting to be a slowly- varying func- 
tion of redshift, such as for example, a downweighting 
variations at z = 4 relative to those at z — 1. It may be 
worth exploring such an optimization in the future. 

Fig. [7] shows d a ii, the mean-subtracted dN/dz for the 
case of all 20 SED templates. The blue (solid) line is 
the actual residual and the cyan (dashed) line is the re- 
construction from the first three principal components. 
While some of the detailed structure is missed, the over- 
all shape is reproduced. The addition of more modes 
would continue to improve the reconstructed d a ii- 

What about the variance of the variation? In the top 
panel of Fig. [8] we plot the standard deviation of the set 
of d a at each point in redshift space. This essentially 
captures the uncertainty in the photometric redshift cal- 
culations as a function of redshift. However, it is not 
necessary to retain all the modes to capture this level 
of error; rather with just three modes over 80% of this 
uncertainty can be represented. This is illustrated in the 
bottom panel of Fig. [8] where we reconstruct the d a only 
using modes 4 and higher and recalculate the standard 
deviation as a function of z using these reconstructed 
d a . It can be seen that the standard deviation from the 
reconstruction is greatly reduced for all collections - in 
other words, a large portion of the variance is captured 
by the first three principal components and little is ac- 
counted for by the higher order modes. 

4. EXAMPLE 2: MAGNITUDE AND TYPE PRIOR 
DISTRIBUTION UNCERTAINTY 

Another modeling assumption that can affect the pho- 
tometric redshift estimations is the assumed distribution 
for redshifts as a function of type and magnitude. In 
the Bayesian framework, these distributions can be ap- 
plied as a prior in the SED template fitting method with 
the intent to i mprove photometric redshift estimates. In 
(|Benitezll2000l) these priors are derived semi-empirically, 
that is to say, a functional form is chosen and the param- 
eters are fit with a training set of spectroscopic redshifts. 
The probability of a galaxy being at a certain redshift is 
given by 

p{z\C,m) <x^2p{z,T\mo)p{C\z,T) (23) 

T 

and using the product rule 

p(z,T\m a ) =p(T\m )p(z\T,m a ). (24) 

The distributions, p(T\mo) and p(z\T, mo) are con- 
strained using auxiliary data. Errors may be introduced 
into these priors in a number of ways. For instance, a 
biased spectroscopic redshift sample (one that underrep- 
resents a certain galaxy population present in the pho- 
tometric sample) could cause errors in the inferred type 



distributions. An incorrect analytical form could also be 
chosen that is unable to accurately represent the prior 
distributions; or simple statistical errors, such as those 
resulting from the finite training set, could cause the 
prior parameters to be known only to low accuracy. In 
this section, we show again that PCA can be used to de- 
scribe uncertainties in redshift estimates resulting from 
such errors in the magnitude and type priors. 

4.1. Mock Data 

To investigate the effects of imperfect knowledge of 
the Bayesian prior, we generate a mock "training" set 
of 100, 000 galaxies. The simulated data are created ac- 
cording to the procedure described in §|3j the generat- 
ing r-band magnitudes, types and redshifts via Monte 
Carlo, while following the distributions for p(T\mo) and 
p{z\T, mo), given in Equations 1161 and 1171 with the free 
parameters of the distribution given in Table [TJ 

Instead of using the LRT template set as a ba- 
sis for generating the 17 band photometry, we use 
the CWW+KIN set that are included with the BPZ 
software, assigning elliptical/reddest as type 1 and 
SB2/blucst as type 6. Rather than limiting ourselves to 
an overly simplified model with only six discrete types, 
which would result in very few type misidentifications, we 
create a continuous color distribution by linearly inter- 
polating between the galaxy types (a continuous number 
between 0.8 and 6.2) to assign color to each galaxy. 

4.2. Results 

We use the mock catalog as a starting point for probing 
the uncertainty in the prior parameters, ft, k t , Zo t and 
at, and investigating the effects of this uncertainty on the 
photometric redshift estimations. We run BPZ to fit each 
galaxy to one of the standard 6 HDFN templates. (Be- 
cause of the interpolation step (described above), some 
galaxies will necessarily be mistyped and these errors, 
on top of statistical errors associated with a finite train- 
ing set, could further affect the estimation of the prior 
parameters.) 

With this output, we use the Monte Carlo Markov 
Chain (MCMC) method (Metropolis Hastings algorithm) 
to find the best fit values and the confidence bounds on 
/{, k t , ZQ t and at- (Note that although there are many es- 
timators that may have worked in our case, MCMC was 
chosen because of the high dimensionality of the parame- 
ter space: 14 free parameters.) The results are presented 
in Table [2] It can be seen that these values are close to, 
but not the same as, the fiducial parameters in Table [TJ 

TABLE 2 

Parameter values for fit using MCMC for magnitude and 

TYPE PRIORS 



Spectral type 


ft 


at 


zot 




kt 


E/SO 

Sbc,Scd 

Irr 


0.305 
0.503 
0.192 


2.350 
1.640 
0.902 


0.416 
0.356 
0.0687 


0.0935 

0.070 

0.119 


0.368 
0.152 



To probe the effects of variations in the priors, we select 
100 values of the parameters drawn from the Markov 
chain. To make sure we consider both best and worst case 
scenarios, we also add the best fit values and parameter 
values at the upper and lower 2a confidence levels (note 
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that when all the parameters are at the 2a confidence 
limit, this reflects a much larger than 2a fluctuation.) 

We run the photo- z redshift code, BPZ, on our simu- 
lated galaxy catalog for each instance of prior parame- 
ters. As in the template uncertainty case, for each galaxy 
in our catalog, the product of the code is a probability 
distribution, P(z), for < z < 4. 

As for the LRT template variations, we sort the galax- 
ies into collections by defining 5 redshift ranges evenly 
spaced from z = to z = 4 and by using the criteria 
that the peak of the P(z) fall within one of these ranges 
to assign each galaxy to a given collection or "redshift 
bin" . As for the previous case, the dN/ dzi for each col- 
lection, i, is calculated by summing the P(z) for each 
galaxy in the collection. 

Variations in the priors create changes in the esti- 
mated dN/dz. The dN/dzi and the residuals, dN/dzi — 
c?7V/dzbost-fit is plotted in Fig. |H1 The differences between 
the overall dN/dz's plotted here and those of the LRT 
model (shown in Fig. [4]), are due to the mock data and 
the template set used for the photo-2 estimation (not 
the priors). While we avoid the overly simplistic 6 tem- 
plate (CWW+SB) model by generating data to include 
linear combinations of the SEDs, we only compare to 
the six discrete HDFN templates to estimate the red- 
shifts. Thus, not every galaxy is accurately represented 
in the template library. It is the non-representativeness 
of these templates that cause the deviations seen in 
Fig. 13 (Note that in the LRT case, we left one tem- 
plate out of 20 instead of using 6 to represent a contin- 
uous mix of types.) The variations in the dN/dz due 
to variations in the prior are small, as can be seen by 
the difference between the dN/ dz using the best fit prior 
parameters and the dN/dz's from each Monte Carlo sam- 
pling dN/dzi — dN/ dzhest-Bt- They are a factor of 5 or 
more smaller when compared to the variations due to 
the template uncertainty. Thus, at least for the case of 
a well-sampled training set, small errors in the prior are 
not as important as having well-calibrated SED template 
set. 

Once again, we perform PCA on the set of dN/dz to de- 
rive a lower dimensional parameterization of the photo- 
z uncertainties. We find that for variations of the priors, 
only two modes are necessary to account for over 80% of 
the variance. The top panels of Fig. [TOlshow the first two 
principal components for each collection. The bottom 
panels show the eigenvalue spectrum for each collection. 

One caveat to this approach is that we have not in- 
cluded the differences in the dN/dz's resulting from the 
difference between the fiducial prior parameters and the 
best fit parameters. However these differences in the 
dN/dz remain small and are of the same order as the 
variations we do consider. We note that in a case where 
the type distributions of the training set are more dras- 
tically biased, the prior parameter distribution may have 
a greater impact. What we have presented here is a 
method for exploring such uncertainties and parameter- 
izing their effects on the redshift distributions in such 
a way that they could be propagated through to error 
bounds on cosmological parameters. 

5. CONNECTION TO OTHER WORK 

Many authors have considered the propaga- 
tion of photo-z errors through to the cosmologi- 
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Fig. 9. — Continuous HDFN priors: Variations in prior param- 
eters drawn from Markov chain. The top panels show the dN/dz 
per collection (tomographic redshift bin), the bottom panels show 
the difference between the dN/dz from the best fit values and the 
dN/dz from the variations, dN/dz^f — dN/dzi, for each collection. 
These differences are small in comparison to those of the LRT tem- 
plate variations. 
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Fig. 10. — Continuous HDFN priors: Top panels of the plot show 
the first two principal components for each collection. The bottom 
panels show the eigenvalue spectrum (i.e. amount of variance cap- 
tured by each principal component). In all cases, more than 80% 
of the variance is captured by the first 2 principal components. 



cal p ar ameters, see fo r instance iHuterer fc Takada 
(2005); IMa et all (12001; lAlbrecht et all (120061 ); IZhar] 
(2006 ft: lAbdalla et all f2008D : IBridle fc Kind (|2007D : 
IMa fc Bernstein! (120081). A commonly quoted result 
is from IMa et al.l (|2006t ). where it was shown that for 
a photo-z distribution modeled as a Gaussian with 
redshift dependent mean and standard deviation, that 
for the a two parameter dark energy model, (the 
dark energy density Ode, its equation of state today 
w = Pde//Ode|z=o), the la errors on z hias (z) and a z (z) 
must be less than .01 per "bin" in order to avoid a 
degradation in the dark energy parameters of more than 
50%. Note that this result does not apply to the mean 
and bias for the tomographic redshift bins (collection 
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in the terminology of this paper), rather the analysis 
refers to the mean and bias on "micro bins" of width 
Sz = .1. Moreover, it does not include the powerful 
self-calibration enabled by a joint analysis of WL and 
BAO, due to the shared large scale stucture. 

To put in perspective the variations we produce via 
changes in template sets and magnitude priors, we com- 
pare the size of the resulting variations in the d N/ dz~, per 
collection 7 to those of the model analyzed in iMa et all 
IpOOfih . The dN/dz 1 per collection is given by 



dN 



p h dN 

. , dzp h — (z)p(z ph \z). 

(7) az tot 

P h 



where 



dN 



z a exp[-(z/z f] 



dz tot 

The photo- z distribution is modeled as 



p(z ph \z) = 



1 



2n<j; 



exp 



(z Zph ^bias) 



2 1 



(25) 



(26) 



(27) 



and is defined for "micro" bins of width Sz = .1, each 
with two photo- z parameters the bias, Zbias(-z), an d the 
scatter, a z (z). The dN/dzj for each "macro" bin, or 
collection, is obtained by interpolating the between each 
p(z p h|z)j, to return a continuous function defined for the 
entire redshift range, and then integrating Eq. 1251 

To calculate the fluctuations in each collection result- 
ing from la errors of .01 on Zb ias (z) and cr z (z) for each 
"micro" bin, we start with the photo- z model described 
above. We assume a fiducial case with 40 "micro" bins 
between z = and z — 4 with Sz{z) = .1, Zbias(^) = 
and cr z (z) = .05(1 + z). We then randomly draw 21 sam- 
ples (consistent with the fact that we generated 21 sam- 
ples from variations of template set (we could have used 
103 to correspond to variation of the prior parameters, 
but given that the distribution is Gaussian, the extra 
variations add little information) from a 80-D Gaussian 
distribution with ouas = v scatter = -01. We linearly in- 
terpolate to create a continuous function in z, P{z p \ 1 \z), 
and we perform the integral in equation [55] each time to 
obtain a set of 21 dN/dz a for each collection. 

We define, as in §[3l d a — dN/dz a — (dN/dz), i.e. 
the mean-subtracted data. The dN/dz a and the d a are 
shown in Fig. [TTJ It can be seen that d a are of the same 
order (although somewhat larger) than those of the tem- 
plate variations (Fig. 2]). However, there are no catas- 
trophic "islands" , i.e. non-zero values for the dN/dz far 
outside of the redshift range for which the collection is 
defined. This is also reflected in Fig. [T2] which shows 
the standard deviation in the mean redshift (Eq. \19\ 
and the variance (Eq. [20)) for each collection. What 
is perhaps most interesting to note is the difference be- 
tween the Gaussian case and the model-based variations 
in collections 4 and 5, and to a lesser extent collection 1, 
which are subject to the redshift degeneracies associated 
with faint galaxies, as discussed in §[3l In particular, 
these differences are apparent in the standard deviation 
in the variances (bottom panels of Fig. [T2")) . where the 
values are significantly higher for the LRT template set 
variations and even, though to a lesser extent, for the 
variations the HDFN priors. 




Fig. 11. — Gaussian model (from IMa et al.l 120061) with 8a = 
8fi = .01 The top panel shows the dN/dz per collection (to- 
mographic redshift bin), the bottom panels depict the residuals, 
dN/dz f it i — dN/dzi, in each collection. Despite stringent con- 
straints on the 5o~ and 8fi for each "microbin", the residuals in 
the dN/dz resulting from variations at this level are larger than 
those of both model-based variations we consider as shown in Fig. 
[4]and Fig. [9] Catastrophic failures with significant probability well 
outside their expected collection are also absent from this param- 
eterization. 

Given the differences between the Gaussian case and 
our model-based variations, we refrain from drawing any 
conclusions about what this would mean for constraints 
on the dark energy, leaving the propagation of such un- 
certainties through to the cosmological parameters for 
future work. However, it is important to note that we 
have considered a level of uncertainty in the photometric 
redshift errors that is significant for constraints on the 
dark energy. The differences from the Gaussian case un- 
derscore the importance of considering realistic sources of 
uncertainty (although the two variations presented here 
are somewhat simplified examples and not intended to 
be taken as a full exploration of such uncertainties). We 
note that we have parameter ized uncer t aintie s of similar 
order to those considered in IMa et~ l. (2006), but more 
generally and with fewer parameters (15 vs. 40). 

6. DISCUSSION AND CONCLUSIONS 

With the potential to significantly improve our knowl- 
edge of the expansion history of the universe, future ob- 
servational programs are necessary for making theoreti- 
cal progress on the puzzle of dark energy. However sys- 
tematic errors have the potential to wreak havoc on these 
results. In particular, photometric redshift methods, es- 
sential to many future experimental programs, rely on 
modeling assumptions that can cause both biases and 
catastrophic errors in the redshift estimates. In order 
for future observations to deliver on their promise, these 
errors must be reduced or well-understood. 

In this work we have presented a method of explor- 
ing and parameterizing modeling uncertainties associ- 
ated with photometric redshift estimates. We consider 
two sources of modeling uncertainty - the SED template 
set and magnitude priors - and we show that both cause 
errors in the photometrically estimated dN/dz. Though 
both cases utilize simplified models to probe their im- 
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Fig. 12. — Comparison of the standard deviation of the mean and 
the variance of the dN/dzi, the redshift distribution per collection, 
for the three models considered (LRT, priors (Continuous HDFN) 
and Gaussian with priors on mean and spread of .01). It can be 
seen that the variance in the means for the Gaussian model is often 
larger than that of the model-based variations. The larger variance 
in the LRT variations and the prior variations for collections 4 an 
5 is indicative of catastrophic outliers which the Gaussian param- 
eterization cannot represent. 



pact, the results indicate that template selection effects 
are dominant to those of the magnitude priors (mean 
subtracted data, d a — dN/dz a — {dN/dz), are nearly an 
order of magnitude larger for the template variations). 
This result comes with the caveat that the underlying 
parameterization of the prior distribution was perfectly 
known for the template variations and we only consider 
the effects of uncertainties in prior parameters (differ- 
ences could be larger in the case where the functional 
form of the distribution is incorrect or training set is bi- 
ased in some way.) Additionally, we vary the template 
set and the priors independently, although the effects of 
these variations may be correlated. A more sophisticated 
model would enable us to examine such modeling uncer- 
tainties jointly. 

We have shown that in both cases, the variations in the 
dN/dz can be characterized by a principal component 
analysis, which selects out the directions of maximum 
variance. Using the principal components as the new 
basis, the dN/dz can be parameterized by the weights, 
bi, where the dN/dz.i for each collection is constructed 
according to Eq. [9l The uncertainties in the model, 
can then be described in terms of a covariance matrix 
for the bi and can be used in a likelihood analysis for 



weak lensing or other cosmological observations involving 
photometric redshifts. 

PCA allows for dimensionality reduction (here we find 
that two or three modes are sufficient to account for 80% 
of the variation in the data). We note that the precision 
of the reconstruction can be easily "tuned" by retention 
of a larger or smaller number of modes. The exact num- 
ber of modes needed for a given level of fidelity to the 
original data will depend on the exact data set and mod- 
eling uncertainties considered. However it is a simple 
matter to choose a desired level of variation captured by 
the parameterization: X a /^2i\i where Aj is the eigen- 
value associated with a given mode. For the case of five 
redshift bins (or collections) and three principal compo- 
nents, the uncertainty in the photometric redshifts can 
be captured by 15 numbers (a computationally tractable 
result). Moreover, this parameterization is completely 
general and does not resort to ad hoc forms such as Gaus- 
sians. Thus we present a way of parameterizing realistic 
modeling uncertainties and propagating these uncertain- 
ties through to the cosmological parameters. 

More generally, this method of estimating photo- z un- 
certainty and propagating to cosmological parameters 
could be applied to the more realistic observational 
case where two or more probes of geometry and struc- 
ture formation (such as cosmic shear tomography and 
baryon acoustic oscillations) are combined. This cross- 
calibration removes degeneracies, significantly reducing 
sensitivity to systematic error. Finally, there are statis- 
tical calibrations that can be done in regions of the sky 
where there are spectroscopic samples, and this informa- 
tion can be imported to the method we have developed. 
In such cases the mean redshift of the galaxy population 
in every collection may be calibrated by cross-correlating 
that sample with a bright spectroscopic sa mple in angle 
and redshift ((Matthew s fc Newman! [201 Of ). Leveraging 
the new deep-wide spectroscopic surveys, one could go 
beyond a calibration of the mean redshift of that sam- 
ple to knowledge of the distribution P\t{z) of galaxies 
in that collection as a fuction of type. This information 
could then be incorporated into the prior. 
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APPENDIX 
ESTIMATOR BIAS 



Here we show that the simple estimator Equation Q is unbiased if the prior and likelihood are both correct. 
Assuming an underlying redshift distribution dN/dz, we get a posterior distribution for the ith galaxy 



where A(f^ i ) is a normalization factor so that 



J dzPi(z) = 1. 



(Al) 



(A2) 
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For simplicity we have suppressed t & mi, which can be integrated out via Equation ([3]). In what follows, we also drop 
the subscript A and superscript d of / for convenience. We group galaxies that meet certain criteria, e.g., within a 
volume in color space, together and call the collection a "redshift bin" . Redshift is only one application; the collection 
criteria can be more general. In terms of symbolics, we use / G F to denote the criteria. On the one hand, dN/dz in 
this F collection is 

dN F f „dN ,„ s / n \ 

-^ = l f J f ^ P ^ Z) - (A3) 

On the other hand, the estimator is 



dN f 



^ P{z\h). (A4) 



dz ^ 

fi£F 

Calculating the expectation value for this estimator is the same as calculating its value in the limit of very large 
numbers of galaxies, assuming that the galaxies follow the redshift distribution dN/dz. In this large- N limit we can 
replace 

with / df I dz* — P{f\z«), 
JfeF 



E 



fiEF 

so that 



dN F 
dz 



J df J dz,^A(l)P{f\z,)P(f\z)^ 



' A m = ^. (A5) 

f£ F CLZ CLZ 

Since the above derivation does not use any specific form of the likelihood function, and since the collection F is 
generic, the estimator is not biased. In reality, the prior dN/dz and the likelihood P(f\z) are not accurately known, 
and systematic errors in them can still cause biases in the results. 
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